rm(list = ls())

library(RColorBrewer)
library(tidyverse)
library(ggpubr)
library(sf)
library(stargazer)
library(haven)
library(readxl)
library(viridis)
library(data.table)

# Figure 1: Types of Revolving-door Official ####
### load data
firm_person_panel = read_csv("~/Dropbox/revolving_door_jmp/replication_files/firm_person_panel.csv")
rd_official = read_csv("~/Dropbox/revolving_door_jmp/replication_files/rd_official_rep.csv") %>% 
distinct()

# find last job post
rd_last = rd_official %>%
arrange(personid, sentence_seq) %>% 
group_by(personid) %>% 
filter(sentence_seq == max(sentence_seq,na.rm = T)) 

# generate data
data = rd_last %>%
select(personid, type) %>%
group_by(type)%>%
summarise(value = n() ) %>%
mutate(type = ordered(type, levels = c("Party", "Government","Military", "People's Congress&CPPCC","Court&Procuratorate", "Youth League"))) %>% 
arrange(desc(type)) %>%
mutate(prop = value / sum(value) * 100) %>%
mutate(ypos = cumsum(prop)- 0.5*prop)

# pie plot ggplot figure
pie = ggplot(data, aes(x="", y=prop, fill=type)) +
geom_bar(stat="identity", width=1,alpha = .5) +
coord_polar("y", start=0) +
geom_text(aes(y = ypos, label = paste(round(prop),"%")), color = "black", size=4) +
theme_void(15) +
scale_fill_viridis_d(direction = -1) +
theme(legend.title = element_blank(),
legend.position = "right") +
xlab("") +
ylab("")

# data generation for time trend
level_vis = rd_last %>% 
#filter(type == "Government" | type == "Party") %>% 
inner_join(firm_person_panel, by = "personid",multiple = "all") %>% 
group_by(level,year) %>%#
tally() %>%
filter(level != "other")%>% 
ungroup() %>% 
mutate(level = if_else(level == "prefecture", "prefecture", level)) %>% 
filter(!is.na(level))

# add levels to level variable
level_vis$level = ordered(level_vis$level, levels = c("central", "province", "prefecture","county", "township"))

# level time trend ggplot figure
level = ggplot(level_vis,aes(x = year, y = n,fill = level,shape = level, color = level)) +
geom_line(alpha = .5) +
geom_point(alpha = .5) +
theme_classic(15) +
xlab("") +
theme(legend.title = element_blank(),
legend.position = "right") +
guides(color=guide_legend(nrow=2)) +
ylab("Number of Officials") 

ggarrange(pie,level)

ggsave( "/Users/zerenli1992/Dropbox/Apps/Overleaf/subsidies_for_sale_2020/figure/time_trend_level.png",
width = 12,height = 4.5)

# Figure 2: Subsidy Provider and Types####
subsidy_df = read_dta("~/Dropbox/revolving_door_jmp/replication_files/program_replication_jop.dta")
# cross sectional
sector_name <- read_excel("~/Dropbox/revolving_door_jmp/replication_files/sector_name.xlsx") %>% 
mutate(ind_code_res_n = as.character(ind_code_res_n)) # sectoral name

breal_sector = c(100,1000,5000,10000)

# sector variation visualization
sector =subsidy_df %>% 
filter(sectorid != "") %>% 
group_by(sectorid) %>% 
summarise(n = n(),
value = round(sum(subsidy_size,na.rm = T)/1000000,2),
average = round(mean(subsidy_size,na.rm = T)/1000000,2)) %>% 
left_join(sector_name, by = c("sectorid" = "ind_code_res_n"))%>% 
ggplot(aes(x = reorder(`Sector Name`, -average),y = average,fill = average),alpha = .5) +
geom_col() +
theme_classic(15) +
scale_fill_viridis_c(alpha = .5) +
theme(axis.text.x = element_text(angle = 45, hjust = 1,vjust = 1),
legend.position = "none") +
#theme(legend.position = element_blank()) +
# scale_fill_viridis_c(option="B",alpha = .8, trans = "log") +
xlab("") +
ylab("") +
ggtitle("Average Value by Sectors (Million RMB)")

sector

# level of provider visualization
level = subsidy_df %>%
mutate(level = str_to_title(level)) %>% 
group_by(level) %>%
tally() %>%
ungroup() %>%
mutate(ypos = cumsum(n)- 0.5*n,
 level_lab = paste(level, "\n(",n,")",sep ="")) %>%
ggplot( aes(x="", y=n, fill=level,fill_t = ),alpha= .8) +
geom_bar(stat="identity", width=1) +
coord_polar("y", start=0) +
theme_void(15) +
theme(legend.position="none") +
geom_text(aes(y = ypos, label = level_lab), size=4) +
xlab("") +
ylab("") +
scale_fill_viridis_d(alpha = .5) +
ggtitle("Levels of Jurisidictions")

# subsidy type visualization
type =subsidy_df %>% 
group_by(type) %>% 
tally() %>% 
ggplot( aes(x=reorder(type, n), y = n,color = n))+
geom_segment( aes(xend=type, yend=0)) +
geom_point( size=4) +
theme_classic(15) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none") +
xlab("") +
ylab("") +
ggtitle("Types of Subsidies") +
scale_color_viridis_c(alpha = .7)

type_level = ggarrange(type,level,ncol = 2)

ggarrange(type_level, sector, ncol = 1)

ggsave( "/Users/zerenli1992/Dropbox/Apps/Overleaf/subsidies_for_sale_2020/figure/subsidy_summary.png",
width = 13,height = 9)

#Figure 3: Sample Illustration####
theme_map_1 <- function(...) {
theme_classic() +
theme(plot.margin = margin(0, 0, 0, 0, "cm"),
axis.ticks = element_blank(),
line = element_blank(),
legend.position = "bottom",
axis.text = element_blank(),
axis.title = element_blank(),
panel.background = element_blank(),
legend.title = element_blank(),
legend.text = element_text(size = 10),
...
)}

# map theme 
path = "~/Dropbox/revolving_door_jmp/replication_files/CHN_adm2.shp"
pref = st_read(path) %>% 
filter(NAME_1 == "Guangdong") %>% 
mutate(red = if_else(NAME_2 == "Shenzhen","Case Group","Control Group") %>% factor()) 

pref = ggplot() +
geom_sf(data= pref, aes(fill = red)) +
theme_map_1()+
#ylim(26.5, NA) +
scale_fill_manual(values=c("red","white"))
pref
ggsave("/Users/zerenli1992/Dropbox/Apps/Overleaf/subsidies_for_sale_2020/figure/case_control_illustrate.png",
 height = 3.5, width = 7)


# Figure 4: Effect Heterogeneity in Different Subsamples####
library(lfe)
library(forcats)
subsidy_df = read_dta("~/Dropbox/revolving_door_jmp/replication_files/program_replication_jop.dta") 

#sectoral analysis
sector_name <- read_excel("~/Dropbox/revolving_door_jmp/replication_files/sector_name.xlsx") %>% 
mutate(ind_code_res_n = as.character(ind_code_res_n))

firm_num = subsidy_df %>% 
group_by(sectorid) %>% 
tally() %>% 
left_join(sector_name, by = c("sectorid" = "ind_code_res_n"))%>% 
filter(sectorid != "") %>% 
mutate(`Sector Name` = paste0(`Sector Name`," (",n,")", sep = "") ) 

sectorid = unique(subsidy_df$sectorid) %>% sort() 
sectorid = sectorid[-1]
coef = vector()
sector = se = vector()

for (i in 1:length(sectorid)) {
df =subsidy_df[subsidy_df$sectorid == sectorid[i],]
m1 = felm(subsidy_l~ pre_hiring+same_place + rd_offic + log(TotalAsset+1) + log(TotalLiab+1)| stkcd+ year+ jurisdiction | 0|stkcd, df ) %>% summary()
coef[i] = m1$coefficients[1,1] 
se[i]=m1$coefficients[1,2] 
}

subsample = data.frame(coef,se,sectorid) %>% 
left_join(firm_num, by = c("sectorid")) %>% 
mutate(`Sector Name`= if_else(is.na(`Sector Name`), "Others",`Sector Name` ) ) %>% 
filter(`Sector Name` != "Others" ) %>% 
#mutate(`Sector Name` = fct_reorder(`Sector Name`, coef , .fun='median' )) %>% 
mutate(color_sign = if_else(coef >0,1,0)) %>% 
filter(!is.na(coef))

sector = ggplot(subsample)+ 
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_point(aes(y = coef, x = fct_reorder(`Sector Name`, coef))) +
geom_pointrange(aes(y = coef, min=coef - 1.96 * se, max=coef + 1.96* se ,x = `Sector Name`, ) ) +
#geom_errorbar(aes(y = coef, min=coef - 1.64 * se, max=coef + 1.64* se ,x = `Sector Name`, ), lwd=.8, width=0 ) +
# facet_wrap(~type+ variable ,scales = "free") +
coord_flip() +
theme_pubclean() +
theme(axis.text.y=element_text(hjust = 1, size =12),
legend.position = "none") + 
xlab("") +
ylab("Effect of Future Revolving Door")

#administrative level analysis
m1 = felm(subsidy_l ~ pre_hiring +rd_offic+ same_place+ log(TotalAsset+1) + log(TotalLiab+1) |year +stkcd + jurisdiction | 0|jurisdiction, subsidy_df %>% filter(level == "province"))
m2 = felm(subsidy_l ~ pre_hiring +rd_offic+ same_place+ log(TotalAsset+1) + log(TotalLiab+1) | year +stkcd + jurisdiction | 0|jurisdiction, subsidy_df %>% filter(level == "prefecture"))
m3 = felm(subsidy_l ~ pre_hiring + rd_offic + same_place+ log(TotalAsset+1) + log(TotalLiab+1) | year +stkcd + jurisdiction | 0|jurisdiction, subsidy_df %>% filter(level == "county"))

# find errors 
m_list = list(m1,m2,m3)
coef_export = function(x){
x$coefficients[1,1]
}
se_export = function(x){
x$cse[1] %>% as.numeric()
}

spec = c( "Province", "Prefecture", "County")

het_vis = tibble(coef = map(m_list, coef_export) %>% unlist() %>% as.numeric(),
 se = map(m_list, se_export)%>% unlist() %>% as.numeric()) %>%
mutate(model =spec)

#visualization
admin = ggplot(het_vis)+
geom_hline(yintercept = 0, color = "red", linetype = "dashed",size =1) +
geom_pointrange(aes(y = coef, min=coef - 1.96 * se, max=coef + 1.96* se ,x = model), lwd=1) +
theme_pubclean() +
theme(legend.title = element_blank(),
legend.position = "none",
legend.text = element_text(size = 12),
axis.ticks = element_blank(),
) +
xlab("") +
ylab("Effect of Future Revolving Door") +
ggtitle("Level of Jurisdiction")


ggarrange(admin,sector)

ggsave( "~/Dropbox/Apps/Overleaf/subsidies_for_sale_2020/figure/het.png",width = 13,height = 5)


